#ifndef __CTCR2D__
#define __CTCR2D__

#include <stdbool.h>
#include "enumerations.h"
#include "constants.h"
#include <stdio.h>
#include <stdlib.h>
// ***********************************************************************
// CR element, conforming, 2D
// ***********************************************************************
static int C_T_CR_2D_dof[3] = {0, 1, 0};
// double C_T_CR_2D_nodal_points[6] = {0.0, 0.0,1.0,0.0,0.0,1.0};
static int C_T_CR_2D_Num_Bas = 3;
static int C_T_CR_2D_Value_Dim = 1;
static int C_T_CR_2D_Inter_Dim = 1;
static int C_T_CR_2D_Polydeg = 1;
static bool C_T_CR_2D_IsOnlyDependOnRefCoord = 0;
static int C_T_CR_2D_Accuracy = 1;
static MAPTYPE C_T_CR_2D_Maptype = CR;

static void C_T_CR_2D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
  // printf("Come to the base D00\n");
  double xi = RefCoord[0];
  double eta = RefCoord[1];
  // printf("get the ref coord!\n");
  values[0] = 2*(-0.5 + xi + eta);
  values[1] = 1.0 - 2.0 * xi;
  values[2] = 1.0 - 2.0 * eta;
  // printf("use the CR interfunction!\n");
  // printf("values = [%2.8f,  %2.8f,  %2.8f], End of base values!\n",values[0], values[1], values[2]);
}

// base function values
static void C_T_CR_2D_D00(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  // printf("Come to the base D00\n");
  double xi = RefCoord[0];
  double eta = RefCoord[1];
  // printf("get the ref coord!\n");
  values[0] = 2*(-0.5 + xi + eta);
  values[1] = 1.0 - 2.0 * xi;
  values[2] = 1.0 - 2.0 * eta;
  // printf("values = [%2.8f,  %2.8f,  %2.8f], End of base values!\n",values[0], values[1], values[2]);
}

// values of the derivatives in xi direction
static void C_T_CR_2D_D10(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  values[0] = 2;
  values[1] = -2;
  values[2] = 0;
}

// values of the derivatives in eta direction
static void C_T_CR_2D_D01(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  values[0] = 2;
  values[1] = 0;
  values[2] = -2;
}
// values of the derivatives in xi-xi  direction
static void C_T_CR_2D_D20(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  values[0] = 0;
  values[1] = 0;
  values[2] = 0;
}
// values of the derivatives in xi-eta direction
static void C_T_CR_2D_D11(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  values[0] = 0;
  values[1] = 0;
  values[2] = 0;
}
// values of the derivatives in eta-eta direction
static void C_T_CR_2D_D02(ELEMENT *Elem, double Coord[], double RefCoord[], double *values)
{
  values[0] = 0;
  values[1] = 0;
  values[2] = 0;
}

static void C_T_CR_2D_Nodal(ELEMENT *elem, FUNCTIONVEC *fun, INT dim, DOUBLE *values)
{
  // CR元的NodalFun对应的是函数在每条边上对应的积分 下面计算积分
  // 首先生成积分格式
  // QUADRATURE *Quadrature = BuildQuadrature(QuadLine4);

  // 直接给出积分格式
  INT NumPoints = 4;
  DOUBLE QuadX14[4] = {0.069431844202973713731097404888715,
                       0.33000947820757187134432797392947,
                       0.66999052179242812865567202607053,
                       0.93056815579702628626890259511129};
  DOUBLE QuadW14[4] = {0.17392742256872692485636378,
                       0.3260725774312730473880606,
                       0.3260725774312730473880606,
                       0.17392742256872692485636378};

  double coord[2];
  double *quadvalue = malloc(dim * sizeof(DOUBLE));
  int i, j, k, m, n;
  // 将边的积分值初始化为0
  memset(values, 0.0, 3 * dim * sizeof(DOUBLE));
  // 对自由度进行循环
  for (i = 0; i < 3; i++)
  {
    m = (i + 1) % 3;
    n = (i + 2) % 3;
    // 对积分点进行循环
    //  for (j = 0; j < Quadrature->NumPoints; j++)
    for (j = 0; j < NumPoints; j++)
    {
      // 计算积分点坐标
      //  coord[0] = Quadrature->QuadX[j] * elem->Vert_X[m] + (1 - Quadrature->QuadX[j]) * elem->Vert_X[n];
      //  coord[1] = Quadrature->QuadX[j] * elem->Vert_Y[m] + (1 - Quadrature->QuadX[j]) * elem->Vert_Y[n];
      coord[0] = QuadX14[j] * elem->Vert_X[m] + (1 - QuadX14[j]) * elem->Vert_X[n];
      coord[1] = QuadX14[j] * elem->Vert_Y[m] + (1 - QuadX14[j]) * elem->Vert_Y[n];
      // 计算积分点对应函数值 累加到边的积分上
      fun(coord, dim, quadvalue);
      // values[i * dim] += quadvalue * Quadrature->QuadW[j];
      for (k = 0; k < dim; k++)
        values[i * dim + k] += quadvalue[k] * QuadW14[j];
    }
  }
  free(quadvalue);
}
#endif
